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A numerical algorithm for studying strongly correlated electron systems is proposed. The 
groundstate wavefunction is projected out after numerical renormalization procedure in the 
path integral formalism. The wavefunction is expressed from the optimized linear combination 
of retained states in the truncated Hilbert space with numerically chosen basis. This algorithm 
does not suffer from the negative sign problem and can be applied to any type of Hamiltonian 
in any dimension. The efficiency is tested in examples of the Hubbard model where the basis of 
Slater determinants is numerically optimized. We show results on fast convergence and accuracy 
achieved with small number of retained states. 
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Numerical methods for studying interacting quantum 
particles have been explored extensively in these decades 
under the rapid progress of computer power. In partic- 
ular, correlated electron systems have been the subject 
of continuing intensive studies. Useful and efficient al- 
gorithms for interacting fermions contribute to under- 
standing of many fundamental issues in condensed mat- 
ter physics as well as to getting insight and prescriptions 
for future potential applications of electronics. However, 
since the proposal of FeynmanEP , simulation methods of 
interacting fermions have remained challenging due to 
several reasons. 

Quantum Monte Carlo (QMC) method is one of the 
powerful techniques for correlated quantum systems and 
it has indeed made possible to clarify various aspects 
of strong electron i^pixelations near the Mott insulator 
in two dimensionsB aai 1 . However, for further efficient 
simulations, the QMC method on fermion sjcstems is 
known to suffer from the negative sign problemcP . If the 
sign problem arises after sampling or truncation of the 
full Hilbert space, it leads to exponentially damped sig- 
nals for increased system sizes and lowered temperatures 
under the cancellation of positive and negative terms, 
thereby makes reliable estimate practically impossible in 
the presence of statistical or round-off errors. 

Density matrix renormalization group developed by 
White© is another technique for fermion systems. It 
offers an efficient way to reach the ground state wave- 
function as well as thermodynamic quantities without 
the sign problem. However, this method can be applied 
efficiently only to systems with one-dimensional configu- 
rations under open boundary conditions. 

In this letter, we propose a new alternative algorithm 
to compute the ground state properties. In our method, 
the optimized ground-state wavefunction |$) is obtained 
as a linear combination of states as 1$) = X^z^lvz) 
within the allowed dimension, L, of the Hilbert space in 
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a numerically chosen basis {| The ground state is 

projected out after successive renormalization processes 
in the path integral. Our renormalization method opti- 
mizes both the basis \ipi) and the coefficients cj. From its 
formalism, it is apparent that this method is completely 
free from the sign problem, because the explicit form of 
the ground state wavefunction is constructed. Physical 
quantity A can be calculated from (A) = 
This may be viewed as a numerical procedure to find the 
best variational ground-state wavefunction within the al- 
lowed dimension of the Hilbert space, L. With increasing 
L, |\&) can be systematically improved from the starting 
variational state at L = 1 chosen as such as the Hartree- 
Fock and RVB states. 

In principle, the renormalization is achieved by succes- 
sively operating the projection operator exp[— tH] with a 
finite r or (H + c) with an appropriate constant c to the 
initial trial wavefunction |$o). The operations project 
out higher energy states and the lower and lower energy 
states are obtained. It is viewed as the numerical renor- 
malization in the imaginary time direction. After operat- 
ing exp[— tH] or (H + c), the dimension of the obtained 
state in our chosen basis of the Hilbert space expands 
through the off-diagonal element of H in the given basis 
representation. In the process of successive operations, 
the dimension exponentially increases. If the system size 
is small enough, the whole Hilbert space can be stored in 
the computer memory. However, with increasing system 
size, of course, the whole space cannot be stored and we 
seek for the best truncation of the Hilbert space within 
the allowed memory and computation time. Under the 
constraint that the number of the stored states, L, is 
kept constant, we repeat the process of the projection 
and the truncation to lower the energy of the resultant 
wavefunction. 

For the renormalization, we repeatedly operate 
exp[— tH] with small r or (H + c) to obtain the ground 
state as |$) = exp[-ri?]P|$ ) or |$) = (H + c) p |$ ) 
for large p. In this paper, we take the first choice 
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cxp[— tH] for the projection. The iterative renormaliza- 
tion process is performed in the following steps. When 
= J2i=i c i Ifi^} is given from the previous (p—1)- 

st step, |* p ) = Et=iEU c ? ) W!? 1) ) in thc P~ th 

step is computed where the set ipj p is provided from 

Ej=i lW = cx p[~ T H] W'f' 1 ) an d J denotes each 

term in the summation over the space expanded through 
the operation exp[— tH] in the low order of r. Thus the 
dimension of thc Hilbert space is expanded due to the 
summation over j. To keep the dimension at the pth 
step, we next pick up L states, out of the ex- 

panded states ipj P j +1 ^ ■ To select L states out of LJ states 
in -0, we solve a generalized eigenvalue problem 

^ ^ Hm.nCn ^ ^ ^ Fm,n^ni (1) 
n n 

where H m>n = (ip { £ +1) \H\ipl? +1) ) and F mn = 
( t Pm +1 ^\'Pn +1 ^}- Note that the basis functions \<fii) gen- 
erated after the operation of cxp \—tH] is not necessarily 
orthogonalizcd. If so, the eigenvalue problem has to be 
generalized accordingly due to nonzero F m , n for m ^ n. 
For each candidate of thc truncated set {|<£m +1 '')} ; we 
calculate the lowest eigenvalue Ao and compare them. 
The set |<^[ p+1 ^) (I = l...,L) are employed when it gives 
thc lowest Ao among the candidates, \tp\ p ^ rV> ). The coef- 
ficients Cn +1 ^ are automatically given from the eigenvec- 
tor with the lowest eigenvalue Ao in the above generalized 
eigenvalue problem (1). Then the (p+ l)-st renormalizcd 
state is given as = c[ p+1) \ip\ p+1) } . 

There may exist several ways to pick up L states out 
of LJ states. When the whole exp[— tH] is operated, 
the expansion J in general becomes huge for large sys- 
tem sizes, because each local term of the Hamiltonian 
may generate off-diagonal terms. Then the number of 
truncating ways becomes exponentially large. To avoid 
it, several ways are possible. One possible way is to 
rewrite exp[— tH] as ]J Li cxp[— rHi] with local terms Hi, 
and perform the truncation process after each operation 
of the local term exp[— rHi]. We call this the local al- 
gorithm. Another way is to introduce a stochastic sam- 
pling, where, among huge number of states expanded 
after the operation of exp[— tH], sampling and anneal- 
ing are performed using the Monte Carlo method and 
one sample is employed after annealing. We call it thc 
Monte Carlo truncation algorithm. 

To present our algorithm in a more concrete way, we 
take the Hubbard model defined by 

H = H t +J2 H m- (2a) 
H t = -tJ2( c l^ + h.c.) (2b) 

(ij) 

and 

Hm = U(n^-^)( nil -^), (2c) 

where M = ^ ia n ia . Here, we follow the standard 
notation. In our example of the Hubbard model, we 
employ, as the chosen basis, Slater determinants, \ip) 



whose matrix representation has N by M matrix for M 
fermions in an TV-site system. We first prepare linearly 
independent L Slater determinants (I = l,...L) 

and their coefficients cf^ to start with the initial state 
|$o) = J2i ^l^i^)- To realize faster convergenccit is 
useful to take, as the initial trial Iv 5 ^), for example, a 
set of samples generated in the conventional quantum 
Monte Carlo calculations. 

We obtain the optimized wavefunction, |$) = 
J2 l Ci\ifii), after sufficient number of rcnormalization pro- 
cess applied to this initial wavefunction. We note that 
the best Hartree-Fock result should be reproduced at 
L = 1 and increasing L systematically improves it. 

To operate exp[— tH] to a Slater determinant, we first 
take the path integral formalism and exp[— tH] is ap- 
proximated by exp[— r(H t — ^M)] exp[— r J2i Hm] for 
sufficiently small r. Then we use the Stratonovich- 
Hubbard transformation for the interaction part. The 
interaction part exp[— THjji] is replaced with thc sum 
over the Stratonovich variable s as 

exp[-rH m } = ^ V cxp[2as(n, T - n %l ) 
1 s =±i 

Ut 

-— Kt + ^i)] (3) 

where a = tanh -1 tanh(^)']. Because of the sum- 
mation over the Stratonovich variables, the number of 
states to be kept after the operations of exp[— tH] expo- 
nentially increases with increasing number of operations. 

Then truncation of thc states is performed follow- 
ing the above mentioned procedure. When the kinetic 
term exp[— r(H t — ^M)] is operated, the dimension in 
the Hilbert space does not increase. The dimension in- 
creases from L to L + 2 when thc local interaction term 
cxp[— rHui] is operated to |<^| P ^)- Here, we newly add 
two states obtained from the operation of Eq.(3) to 
in addition to thc orig inal |^ (p) ) itself. This increases the 
whole number of states from L to L + 2. Thc original 
\<p^) is also retained as a candidate simply to reach 
a better estimate. In the local algorithm, the trunca- 
tion from L + 2 to L is taken by solving the general- 
ized eigenvalue problem every time at each operation of 
cxp[— rHui] to \(p\ p ^} and by finding the eigenstates with 
the lowest eigenvalue among the sets of L retained states. 
This projection and truncation processes are repeated 
NL times to complete a unit operation of exp[— tH] to 
thc L retained states. In the Monte Carlo truncation al- 
gorithm, JT^ exp[— rHui] is applied to |<^| P ^), which gen- 
erates new 2 N terms. If this projection is performed from 
p = 1 to p = L, it generates new 2 N L states. In this case, 
2 N corresponds to J denned before. For truncation, one 
sample with L states is employed among 2 N L states af- 
ter Monte Carlo sweeps of the Stratonovich variables to 
pick up lower energy states. The convergence to the opti- 
mized states seems to become faster when the above two 
procedures are combined and taken alternately. Each 
stored state \tpi) can be viewed as analog of each sample 
in the QMC method. However, an important difference 
is that the present method is free of the sign problem 
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because the stored L states are used to generate a many- 
body wavefunction by the optimized linear combination 
of them. Another difference is that the coefficients of the 
states, ci, in the linear combination are also optimized for 
a better approximation of the ground state in contrast 
with the QMC method where the weight of each sample 
is taken unity in the importance sampling algorithm. 

The ground state is approximated by a truncation of 
the Hilbert space also in the DMRG method. In the 
DMRG method, the truncated basis is represented after 
many linear transformations of the original basis. Al- 
though this transformation makes it possible to reach 
highly accurate results, the transformed basis is not 
represented easily from a simple basis we usually take. 
This is the reason why the application of DMRG is re- 
stricted to the cases with the one-dimensional configu- 
ration. Otherwise the operator of the physical quanti- 
ties are not easily represented in the transformed ba- 
sis. In contrast with DMRG, our algorithm allows basis 
only when the matrix elements can be easily estimated. 
This constraint would make the truncation error large for 
fixed L, while it allows calculations in any type of lattice 
structures and boundary conditions. Truncations of the 
Hilberi^pacej-SKere also examined in several other pro- 
posalsoH'B WliiP. However, to reach faster convergence 
at smaller L, it is crucial to carefully optimize the wave- 
function iteratively within a fixed L as in the procedure 
of the present work. 

We summarize our renormalization process in the fol- 
lowing. 

(1) dnitial state 

A set of truncated basis \<fi ) (I = 1,L) and the coef- 
ficients Cj are given to generate the initial variational 

wavefunction |$ ) = J2i c i \Pi ) ■ 

(2) : Projection 

\<Pl ) are evolved by the operation of exp[— tH\. After 
this evolution, the number of states to be kept is in- 
creased and thus the dimension in the Hilbert space is 
expanded. 

(3) : Truncation 

From the increased states, we select the same num- 
ber of states, L, as the beginning and discard the rest 
to keep the dimension in the Hilberd space. This selec- 
tion is numericaly performed to realize the lowest energy 
states. In the numerical process, such lowest energy state 
is searched after the numerical optimization of the coeffi- 
cients Ci by examining proposed ways of the truncations. 
In each possible way of the truncations, r, the lowest 
eigenvalue Ao and eigenvector q for a generalized eigen- 
value problem (<E> r |7J |<& r ) = A($ r |$ r ) is solved and Ao's 
are compared to select the lowest energy state. 

The steps (2) and (3) constitute a unit renormalization 
procedure. 

(4) : Iteration 

The chosen sets \<Pt ) and cp (I — l,...,L) are used 
as the starting point of the next-step renormalization 
procedure. 

(5) : Convergence 

This renormalization process is repeated until the en- 
ergy converges to the lowest energy under the allowed 



dimension of the Hilbert space. 

Since we are able to treat only a subspace of the whole 
Hilbert space in large systems, it is important to analyze 
the extrapolation procedure using the calculated results 
under constrained Hilbert space. At the moment, it is 
not clear enough how physical quantities converges to the 
exact value as a function of L. However, it can be shown 
that the quantities linearly converge as a function of the 
energy variance defined by = ({E 2 ) — {E) 2 )/{E) 2 if 
A E is smalM, where (E 2 ) = ($ r |# 2 |$ r )/($ I .|3v) and 
(E) = ($ r \H\$ r )/{$ r \$ r ). Note that the variance van- 
ishes if the exact ground state is obtained. The conver- 
gent result in the ground state may be obtained after 
extrapolation to the zero variance, although we have to 
be careful in excluding possible accidental convergence 
to an excited eigenstate with increasing L. 



Hubbard, 6*2 lattice, 5 up, 5 down electrons 
fully periodic, U=4,t=1 
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Fig. 1. The calculated energy of the Hubbard model at t = 1 and 
U = 4 as a function of the energy variance A^; for 5 up and 5 
down electrons on the 6x2 lattice with the periodic boundary 
condition. 



Here we show a few examples for the Hubbard model 
to show the efficiency of this numerical procedure. Wc 
show results at U/t = 4. The first example is the case of 
5 up and 5 down electrons on the 6 by 2 lattice with the 
full periodic boundary condition. The convergence of the 
ground state energy to the exact value is plotted in Fig.l, 
where the exact value is reproduced in 4 digits after the 
extrapolation — > 0. The obtained energy per site is 
-1.4744 and the exact value is -1.4746. The largest L we 
used in this extrapolation is 400. The momentum distri- 
bution < rik >=< c\ck > and the spin structure factor 
S (q) obtained after the extrapolation are compared with 
the exact values in Table I. 

The second example is the case of 13 up and down 
electrons each on the 6 by 6 lattice with the periodic 
boundary condition. The results obtained from L up 
to 400 again are compared with the QMC results. The 
number of states to be retained seems gradually to in- 
crease with the increasing system size if one wishes to 
reach the same accuracy. However, the size scaling for 
extrapolation to L — > oo appears to equally work. The 
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Table I. Comparison of the equal-time spin correlations S(q) and 
the momentum distribution n(q) between the present results and 
the exact ones for the Hubbard model at U/t = 4 on a 2x 6 
lattice with the periodic boundary condition. 

Equal time spin correlation S(k x ,k y ) 



(kx,ky) 


present result 


exact 


(1,0) 


0.03662 


0.03610 


(2,0) 


0.04125 


0.04075 


(3,0) 


0.04161 


0.04160 


(0,1) 


0.1293 


0.1299 


(1,1) 


0.1302 


0.1306 


(2,1) 


0.1339 


0.1335 


(3,1) 


0.1375 


0.1335 


Jomentum distribution n 


(k x ,k y ) 


(kx,ky) 


present result 


exact 


(0,0) 


0.9695 


0.9681 


(1,0) 


0.9620 


0.9601 


(2,0) 


0.9305 


0.9281 


(3,0) 


0.01676 


0.01850 


(0,1) 


0.06892 


0.06806 


(1,1) 


0.04461 


0.04513 


(2,1) 


0.02672 


0.02787 


(3,1) 


0.02233 


0.02281 



obtained energy per site, -L1751 is compared with the 
QMC result, -1.1756±O.OO20). 

Here, we note several aspects on efficiency of this al- 
gorithm clarified in our study and also discuss problems 
to be solved in the future. In the present algorithm, 
solving the generalized eigenvalue problem (1) requires 
the computation time proportional to L 3 . This is the 
time needed to update one of \<p\ p) ) to Then 
the total computation time is proportional to L 4 for up- 
dating all the \tp\ p) ) {I = 1,..,L) once. The evolution 
cxp[— tH] I $( p )) requires LN 2 operations. To compute 
all the matrix elements H m ^ n and F m ^ n , we need time 
proportional to L 2 N 3 because computation of each el- 
ement H m ^ n or F m . n takes time oc iV 3 . Since the ma- 
trix elements can be calculated from the Green func- 
tion Gi } i'(k,m) = ((/S/lcfccJJ^J), it is still proportional to 
L 2 N 3 in updating Y\ i exp[— tHui] \<pi) by the local algo- 
rithm, because exp[— rf/^jjc/ccj^ exp[— rHmW^) can 
be calculated from the old Green function Gy/ (k, m) 
using the local updating procedure used in the QMC 
methodEr with the computation time oc TV 2 . When the 
convergence to a required precission is reached in a fixed 
number of renormalization steps, the total computation 
time is dominated by the larger one between L A and 
L 2 N 3 . This computation time may be compared with 
the time N S N 3 needed for N s steps of the QMC calcu- 
lation. If L can be suppressed less than N s to reach the 
same accuracy, the present algorithm may compete or 
works more efficiently than the QMC method. In fact, L 
can be suppressed to several hundreds in our examples 
to reach the accuracy of three digits in energy, which is a 
typical accuracy of QMC calculations with the order of 
10 4 steps. In fact, the statistical error in QMC is scaled 
by 1/ v^Vs • The advantage of this algorithm is that it 
does not suffer from the statistical sampling error as in 



the QMC method. Furthermore, the present algorithm 
is the only efficient way in multi-dimensional systems if 
the sign problem is serious in the QMC calculations. 

We next discuss interesting future subject for more 
improved algorthm. In our study, a combined algorithm 
of Monte Carlo and local ones are tested for the trun- 
cation. To project out into the ground state, a more 
efficient method of projection and truncation could be 
used. In fact, if each ) is evolved separately from 
1 to L, the projection may become slow possibly due to 
the multi- minimum structure of the energy in the Hilbert 
space. To overcome it, it would be useful to implement 
an efficient algorithm employed in random systems with 
many nearly degnerate energy minima structure. 

When L becomes large, the most time-consuming pro- 
cess becomes solving of the generalized eigenvalue prob- 
lem. It would be helpful if the computation time of 
this part would be substantially reduced from L 4 oper- 
ations. Since the matrix operations have instabilities, a 
simple method like the conjugate gradient method does 
not work and a more careful algorithm is needed. 

In summary, a new algorithm for computing physical 
quantities of correlated fermion systems is proposed. The 
ground state wavefunction is obtained as a linear com- 
bination of states within the allowed dimension of the 
Hilbert space in numerically optimized basis. The opti- 
mization is achieved after the numerical renormalization 
procedure in the path integral form. The results are ex- 
trapolated by systematically examining the dependence 
on the dimension of the Hilbert space. Physical quanti- 
ties show good convergence and accuracy in the calcu- 
lation of the Hubbard model. Since this method does 
not suffer from the negative sign problem and can treat 
any type of lattice structure, it opens a possibility for 
efficient simulations which cannot be reached by existing 
algorithm such as the quantum Monte Carlo method or 
the density matrix renormalization group. 
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